//      arma_helper.h
//      
//      Copyright 2012 tqlong <tqlong@espada>
//      
//      This program is free software; you can redistribute it and/or modify
//      it under the terms of the GNU General Public License as published by
//      the Free Software Foundation; either version 2 of the License, or
//      (at your option) any later version.
//      
//      This program is distributed in the hope that it will be useful,
//      but WITHOUT ANY WARRANTY; without even the implied warranty of
//      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//      GNU General Public License for more details.
//      
//      You should have received a copy of the GNU General Public License
//      along with this program; if not, write to the Free Software
//      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
//      MA 02110-1301, USA.

#ifndef ARMA_HELPER_ONLY_H
#define ARMA_HELPER_ONLY_H

#include <armadillo>

#define vec_t arma::vec
#define matrix_t arma::mat
#define row(m,i) ((m).row(i))
#define col(m,i) ((m).col(i))

struct line {
  int i_;
  int j_;
  double s_;
  
  line() {}
  line(int i, int j, double s) : i_(i), j_(j), s_(s) {}      
  static void write_block(std::ostream& o, line l) {
    o.write((char*)&l, sizeof(line));
  }
  std::istream& read(std::istream& i) {
    return i.read((char*)this, sizeof(line));
  }
};

template<typename MATRIX>
void read_sparse_matrix(const std::string& inFile, MATRIX& A, bool is_binary) 
{
  FileAccessRead(in, inFile);
  if (!is_binary) {
    int m, n;
    double check_sum;
    in >> m >> n >> check_sum;
    A.set_size(m, n);
    double sum = 0;
    for (int i = 0; i < m; i++)
      for (int j = 0; j < n; j++) {
	double s;
	in >> s;
	if ( fabs(s) > 1e-16 ) A(i,j) = s;
	sum += s;
      }
    MSG("check_sum=" << check_sum << " sum=" << sum << " diff=" << fabs(sum-check_sum));
    //if (fabs(sum-check_sum) > 1e-10) FATAL_ERROR("check sum error");
  }
  else { // binary mode
    line head, l;
    head.read(in);
    A.set_size(head.i_, head.j_);
    double sum = 0;
    while (!in.eof()) {
      l.read(in);
      if (in.fail()) break;
      A(l.i_, l.j_) = l.s_;
      sum += l.s_;
    }
    MSG("check_sum=" << head.s_ << " sum=" << sum << " diff=" << fabs(sum-head.s_));
    //if (fabs(sum-head.s_) > 1e-10) FATAL_ERROR("check sum error");
  }
  in.close();
}

template<typename MATRIX>
void save_sparse_matrix(const std::string& outFile, const MATRIX& A, bool is_binary) {
  double sum = arma::accu(A);

  FileAccessWrite(out, outFile);
  if (!is_binary) { // text mode
    out << "\t" << A.n_rows << "\t" << A.n_cols << "\t" << sum << ENDL;
    for (uint i = 0; i < A.n_rows; i++) {
      for (uint j = 0; j < A.n_cols; j++) out << "\t" << A(i, j);
      out << ENDL;
    }
  }
  else { // binary mode
    //#define tr1(x, it) for ( typeof((x).begin1()) it = (x).begin1(); it != (x).end1(); it++)
    line::write_block(out, line(A.n_rows, A.n_cols, sum));
    for (uint i = 0; i < A.n_rows; i++)
      for (uint j = 0; j < A.n_cols; j++)
	line::write_block(out, line(i, j, A(i,j)));
  }
  MSG("check sum=" << sum);
  out.close();
}

template<typename MATRIX>
void read_sparse_matrix_subrows(const std::string& inFile, MATRIX& A, bool is_binary, int begin, int end_exclusive) 
{
  FileAccessRead(in, inFile);
  if (!is_binary) {
    int m, n;
    double check_sum;
    in >> m >> n >> check_sum;
    A.zeros(end_exclusive - begin, n);
    double sum = 0;
    for (int i = 0; i < m; i++)
      for (int j = 0; j < n; j++) {
	double s;
	in >> s;
	if (i >= begin && i < end_exclusive)
	  A(i-begin,j) = s;
	sum += s;
      }
    MSG("check_sum=" << check_sum << " sum=" << sum << " diff=" << fabs(sum-check_sum));
    //if (fabs(sum-check_sum) > 1e-10) FATAL_ERROR("check sum error");
  }
  else { // binary mode
    line head, l;
    head.read(in);
    A.zeros(end_exclusive - begin, head.j_);
    double sum = 0;
    while (!in.eof()) {
      l.read(in);
      if (in.fail()) break;
      if (l.i_ >= begin && l.i_ < end_exclusive)
	A(l.i_-begin, l.j_) = l.s_;
      sum += l.s_;
    }
    MSG("check_sum=" << head.s_ << " sum=" << sum << " diff=" << fabs(sum-head.s_));
    //if (fabs(sum-head.s_) > 1e-10) FATAL_ERROR("check sum error");
  }
  in.close();
}

#endif
